Predicting C-H/tt interactions with nonlocal density functional theory 
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We examine the performance of a recently developed nonlocal density functional in predicting 
a model noncovalent interaction, the weak bond between an aromatic vr system and an aliphatic 
C-H group. The new functional is a significant improvement over traditional density function- 
als, providing results which compare favorably to high-level quantum-chemistry techniques but at 
considerably lower computational cost. Interaction energies in several model C-H/vr systems are in 
generally good agreement with coupled-cluster calculations, though equilibrium distances are consis- 
tently overpredicted when using the revPBE functional for exchange. The new functional correctly 
predicts changes in energy upon addition of halogen substituents. 



I. INTRODUCTION 

Accurate treatment of van der Waals bonding repre- 
sents one of the primary challenges for modern quantum 
chemical calculations. Such interactions are important 
in a great diversity of chemical and biological systems, 
and have been difficult to model in a highly accurate 
and computationally tractable manner. State-of-the-art 
techniques such as coupled cluster theory (CCSD(T)) 
generally treat such interactions very accurately, but 
are computationally demanding and often depend sen- 
sitively on the basis set used. Density functional the- 
ory (DFT) methods are attractive in terms of computa- 
tional efficiency, but common exchange-correlation func- 
tionals used do not include effects such as dispersion in- 
teractions. For this reason, traditional DFT calculations 
frequently provide a poor description of van der Waals 
bonding in molecular systems. 

Here we consider a recently developed van der Waals 
density functional (vdW-DF) which explicitly includes 
nonlocal electron correlation without the use of any em- 
pirical parameters [H |2] . This method has been success- 
fully applied to a number of aromatic tt systems such 
as dimers of benzene [SIS], napthalene, anthracene, and 
pyrene [5], as well as a variety of bulk systems such as 
crystalline polyethylene [6 . The vdW-DF functional was 
strongly influenced by the physics included in its forerun- 
ner [3IH1I1]. 

Of key importance is the ability of this nonlocal func- 
tional to treat a range of different noncovalent interac- 
tions. Along these lines we consider the attractive bond- 
ing between an aliphatic C-H group and an aromatic tt 
system, a configuration which has received considerable 
attention in recent years [inilll]. This so-called C-H/tt 
interaction is frequently referred to as a weak hydrogen 
bond, though recent studies have demonstrated that it 
arises primarily from dispersion, with electrostatic at- 
traction playing a minor role [T^ [T^ [H] . The C-H/tt 
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FIG. 1; The orientation of (a) CHXs-benzene molecules, 
where X=II, F, or CI, and (b) the methane-indole system. 



interaction is important in molecular packing [15j, pro- 
tein interactions [HI [T^, and in a variety of biological 
systems [TSJ [H]. The wide interest in this interaction, 
combined with the availability of high level CCSD(T) 
and MP2 calculations, make it a valuable test case for 
the new van der Waals functional. 

In this paper we consider a number of model sys- 
tems which allow careful study of the C-H/tt interaction. 
These include the methane-benzene system, trifluoro- 
and trichloromethane with benzene, and the methane- 
indole complex. All results are compared with high- 
level CCSD(T) and MP2 calculations. The interaction 
energies from vdW-DF are a significant improvement 
over gradient-corrected and hybrid DFT functionals, and 
compare favorably to CCSD(T) results extrapolated to 
the limit of a complete basis set. Equilibrium distances 
are generally overpredicted by vdW-DF . Our results sug- 
gest that this functional, which requires little time be- 
yond a standard DFT calculation, holds great promise 
for modeling large organic and biological systems. 



II. NONLOCAL DENSITY FUNCTIONAL 

Our approach to calculating the C-H/tt interaction 
uses a recent density functional which explicitly includes 
van der Waals correlations [1,, 2 , 7, 8, 9 . In this frame- 
work, which we will abbreviate as vdW-DF, the total 
energy functional is written as 
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E[p] = [p] + [p] + Eh [p] + [p] + E^ [p] + i?f ^ [p] , 

(1) 

where Tg is the independent-particle kinetic energy, Vpp is 
the ionic pseudopotential, and Eh is the Hartree energy. 
For the exchange-correlation functional, we separate the 
exchange energy E^, the local correlation E^ , and the 
nonlocal component E^^ . Following previous work, we 
use the revPBE expression for exchange [20^ for all calcu- 
lations. E^ is approximated here by the LDA correlation. 
The nonlocal term, which is key for capturing long-range 
correlations such as dispersion, can be expressed in a gen- 
eral form as 



^= = 2 yy *idr2p(ri)0(ri,r2)p(r2). (2) 

The kernel (/)(ri, depends only on the distance |ri — 
r2| and the density p near ri and V2- Details regarding 
explicit evaluation of the kernel are given in Ref. [T]. No 
empirical parameters are used in the nonlocal correlation. 
Though the revPBE exchange does not satisfy the Lieb- 
Oxford bound for any arbitrary density as does PBE, 
in practice it satifies the true, integrated Lieb-Oxford 
condition for typical systems [^Ul [H]. Thus our overall 
approach, using vdW-DF and revPBE exchange, may be 
considered a truly first-principles method for treating van 
der Waals systems. 

In previous work E^^ was computed as a post- 
processing step using the charge density from a normal 
DFT calculation [5J 2] • Explicit evaluation of the non- 
local potential V^^ in a self-consistent way is possible 
through recently derived analytical expressions [2] . How- 
ever, test calculations for the methane-benzene system 
indicated little difference in the interaction energies be- 
tween these two approaches. Thus for all systems we 
continue to utilize the post-processing method, which re- 
quires very little time beyond a standard DFT calcula- 
tion. 

All vdW-DF calculations were done in a modified ver- 
sion of abinit [25]. Norm-conserving TrouUier-Martins 
pseudopotentials were used, with a kinetic energy cutoff 
of 50 Rydberg. Calculations were performed in a large 
cubic cell with at least 18 A of space between adjacent 
periodic images. The nonlocal energy was calculated on 
a real-space grid identical to the one used for Fourier 
transforms in the planewave code. Calculations for hy- 
brid DFT functionals were performed using CRYSTAL03 
i23j. 

Molecular geometries were optimized using the B3LYP 
hybrid functional with an aug-cc-pVTZ basis set; this 
optimization gave structures that were generally within 
0.01 A of previous high-level calculations for methane- 
benzene and methane-indole |14j . These same geome- 
tries were also converged with respect to forces derived 
from the self-consistent implementation of vdW-DF, in- 
dicating that the nonlocal functional vanishes appropri- 
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FIG. 2: The interaction energy of the methane-benzene dimer 
using vdW-DF. Comparison is made with the PBE [27] and 
revPBE [20j GGA functionals, the B3LYP hybrid functional 
[28] . as well as previously reported CCSD(T) and MP2 data 
extrapolated to the limit of a complete basis set (CBS) [14) . 



ately in intramolecular situations. In all cases the iso- 
lated molecule geometry was kept rigid during the two- 
molecule calculations. 



III. RESULTS AND DISCUSSION 
A. Methane-benzene system 

The methane-benzene complex represents one of the 
simplest interacting C-H/tt systems. Following the re- 
sults of previous experimental [53] [^^ and theoretical 
[Ml |26] work, the molecules are oriented such that a sin- 
gle C-H group on the methane points directly at the cen- 
ter of the benzene, as shown in Fig. [T] The separation is 
then defined as the distance between the central C atom 
in methane and the center of the benzene ring. The sepa- 
ration axis is thus coincident with the Cg symmetry axis 
of benzene, and the dimer system overall is set in a Cs^ 
symmetry. Similar to the report by Ringer and coau- 
thors [14], the dimer energy is essentially invariant for 
rotations around the dimer separation axis, varying by 
at most lO^"* kcal/mol. 

In Fig. [2] we present our results for the interac- 
tion energy of the methane-benzene dimer. The vdW- 
DF functional is a significant improvement over other 
standard DFT approaches, yielding an interaction en- 
ergy comparable to the CCSD(T) results of Ringer 
et al p3^. Their coupled-cluster interaction energy 
was extrapolated to the limit of a complete basis set 
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(CBS) using Helgaker's method with large augmented 
correlation-consistent triple- and quadruple-C basis sets 
[H] . Shibasaki and coauthors also report similar theoret- 
ical CCSD(T) results PF. Despite its relative simplicity, 
vdW-DF captures the C-H/tt energy quite accurately as 
compared to these high level calculations. The equilib- 
rium dimer separation predicted by vdW-DF is roughly 
8% larger than the CCSD(T) result, which appears to 
be a general trend of vdW-DF when using the revPBE 
functional for the exchange energy [21 |3l IH . 

Three results using traditional DFT methods are also 
shown in Fig. [2] The Perdew-Burke-Ernzerhof (PBE) 
nonempirical generalized gradient (GGA) functional pro- 
duces a binding that is considerably weaker than that 
given by vdW-DF, though its predicted equilibrium dis- 
tance is only slightly higher than the nonlocal functional 
(see Table |l] for numerical values). Zhang and Yang's 
revised PBE functional (revPBE), which is identical to 
PBE except for a semiempirical adjustment of one pa- 
rameter in the exchange portion [20], produces no bind- 
ing. The popular hybrid functional B3LYP, which uses a 
semiempirical mixing of exact exchange and GGA terms 
[28] , also produces no binding for the system. Recent cal- 
culations by Shibasaki and coauthors show similar results 
for a range of GGA and hybrid functionals [2S]. Thus, 
as we would expect for a dispersion-dominated dimer, 
traditional DFT methods compare poorly to high level 
CCSD(T) and MP2 results. 

Shibasaki and coworkers have recently measured the 
experimental interaction energy of gas phase methane- 
benzene clusters by mass analyzed threshold ionization 
spectroscopy f25' . For comparison with the experimental 
value we must subtract the vibrational zero-point energy 
(ZPE) from our predicted C-H/tt interaction energy. The 
predicted change in ZPE for formation of the methane- 
benzene dimer is 0.296 kcal/mol using MP2 with a cc- 
pVTZ basis set and scaling factors to match the ex- 
perimental and theoretical frequencies for the methane 
symmetric stretching mode [55]. Subtracting this gives 
an interaction energy of —1.27 kcal/mol using vdW-DF 
and —1.15 kcal/mol using CCSD(T); both of these com- 
pare favorably to the experimental binding energy, which 
was determined to be between —1.03 and —1.13 kcal/mol 

MM- 



B. Trifluoro- and trichloromethane 

We next consider the effect of halogen substituents 
on the aliphatic molecule. Some experimental measure- 
ments of the role of such substituent atoms on the C- 
H/tt interaction in solution have been reported, includ- 
ing evidence that benzene forms a stable dimer with 
trichloromethane [301 131] • Here we consider replacement 
of three hydrogen atoms on the methane by chlorine or 
flourine, both of which increase the interaction energy of 
the complex. The results for these dimers of trifluoro- 
and trichloromethane with benzene are given in Fig. [3] 



o 

a 

o 

a 
_o 

03 




vdW-DF 
^ PBE 

• CCSD(T)/CBS_ 

♦ MP2/CBS 



3 3.5 4 4.5 5 5.5 6 
Distance (A) 

FIG. 3: The interaction energy of benzene with tri- 
fluromethane (top panel) and trichloromethane (bottom 
panel). Results are compared to previous CCSD(T) and MP2 
calculations by Tsuzuki and coworkers extrapolated to the 
limit of a complete basis set (CBS) |29) . 



and Table HI 

In both systems the C-H bond points directly at the 
center of the benzene ring, and the remaining three H 
atoms on the methane are replaced. The geometries for 
CHF3 and CHCI3 were optimized independently, and 
the monomer structure was also used for dimer cal- 
culations. Trifluoro- and trichloromethane both have 
increased molecular polarizabilities compared to CH4. 
Methane has a polarizability of roughly 17.9 a.u., com- 
pared to 19.0 a.u. for trifluoromethane and 57.6 a.u. for 
trichloromethane |32j . Since dispersion is the dominant 
component of the C-H/tt interaction, we would expect an 
increased interaction energy for the molecules with larger 
polarizabilities. This is indeed observed, as seen in Fig. [3] 
the benzene-trichloromethane binding is much stronger 
than that seen in the methane system, and is compara- 
ble to the strength of the hydrogen bond in a dimer of 
water. We would also expect an increase in electrostatic 
attraction upon substitution of F or CI, though previ- 
ous results indicate this is a much smaller effect than the 
enhanced dispersion [55^ . 

The performance of vdW-DF in these halogenated 
dimers is similar to that seen in the methane-benzene 
system. The nonlocal functional is a considerable im- 
provement over the prototypical GGA functional PBE, 
providing interaction energies comparable to CCSD(T) 
results extrapolated to the limit of a complete basis set by 
Tsuzuki and coauthors [29]. In these systems the vdW- 
DF interaction energy underestimates the previously re- 
ported CCSD(T) limit, in contrast to a slight overes- 
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TABLE I: Dimer separations (Rd) and interaction energies 
Eint for all C-H/tt systems. Distances are in A and energies 
are in kcal/mol. vdW-DF results are given using the revPBE 
exchange functional. CCSD(T) and MP2 results are taken 
from Refs. [H] and [22]. 
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vdW-DF 


4.1 -1.57 


3.6 -3.73 


3.6 -5.11 


4.0 


-1.81 


PEE 


4.2 -0.403 


3.7 -2.04 


3.7 -2.07 


4.2 


-0.408 


CCSD(T) 


3.8 -1.45 


3.4 -4.20 


3.2 -5.60 


3.8 


-1.57 


MP2 


3.8 -1.79 


3.4 -4.60 


3.2 -7.20 


3.7 


-1.96 



timation for the methane-benzene and methane-indole 
complexes. The hterature values for the latter two sys- 
tems used slightly different basis sets as well as different 
methods for extracting the complete basis limit than were 
used for the halogenated dimers [TU [55] . We note again 
the advantage of using a simple, straightforward DFT 
method such as vdW-DF. 

The vdW-DF functional with revPBE exchange again 
overstimates the equilibrium separation of the dimers by 
approximately 6% for trifluoromethane-benzene and 13% 
for trichloromethane-benzene. MP2 and CCSD(T) cal- 
culations indicate that the latter has a smaller dimer 
separation, as would be expected from the increased 
dispersion component. vdW-DF does not reproduce 
this trend, predicting similar separations for both halo- 
genated dimers. 



C. Methane-indole 

We next consider the interaction of methane with 
a larger aromatic system, the indole molecule. The 
methane is placed over the five-membered ring of indole, 
again with a single C-H bond pointing directly at the cen- 
ter of the ring (see Fig. [ij. Geometries of the monomers 
were optimized separately and kept rigid during dimer 
calculations. In this configuration the interaction energy 
is again essentially invariant to rotations around the axis 
connecting the methane C and the center of the five- 
membered ring. 

Our results for the interaction energy of the methane- 
indole dimer are given in Fig. |4] and Table |lj MP2 and 
CCSD(T) results are taken from the work of Ringer and 
coauthors, where an aug-cc-pVTZ basis set was used 
[14j . The magnitude of the interaction and the dimer 
separation are similar to the methane-benzene dimer, 
though vdW-DF predicts that methane-indole is more 
strongly bound, in agreement with CCSD(T) and MP2 
results. An analysis by symmetry-adapted perturbation 
theory shows a similar distribution of components in the 
methane-benzene and methane-indole dimers, the most 
significant difference being a larger dispersion component 
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FIG. 4: The interaction energy of the methane-indole dimer. 
CCSD(T) and MP2 results are taken from Ref. [H]. 



(0.589 kcal/mol more negative) in the latter. As before, 
vdW-DF overpredicts the CCSD(T) dimer separation by 
approximately 5%. 



D. Conclusion 

We have studied the performance of the van der Waals 
density functional in treating the weak bond between an 
aromatic tt system and an aliphatic C-H group. The 
properties of four prototypical C-H/tt systems were ana- 
lyzed: methane-benzene, trifluoro- and trichloromcthane 
with benzene, and methane-indole. 

The vdW-DF method is a substantial improvement 
over standard DFT functionals, providing in all cases an 
interaction energy that is comparable to high-level quan- 
tum chemistry calculations. In the case of the methane- 
benzene dimer, vdW-DF predicts a total interaction en- 
ergy of —1.57 kcal/mol, compared to —1.454 kcal/mol 
from CCSD(T) and -1.79 kcal/mol from MP2 [14 . Tra- 
ditional DFT functionals predict little or no binding for 
this dimer, and generally wavefunction-based methods 
require large basis sets to achieve values in this range 

The substitution of F and CI for three of the methane 
hydrogens results in an increased binding due to larger 
dispersion and electrostatic interactions. The vdW-DF 
functional correctly reproduces the increased interaction 
energy in good agreement with the CCSD(T) trend, but 
does not match the reduction in equilibrium dimer sep- 
aration. vdW-DF is also able to resolve the small bind- 
ing energy difference between the methane-benzene and 
methane-indole dimers, but overpredicts the equilibrium 
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separation in the latter system as well. 

Overall, the performance of the new van der Waals 
density functional for these very weak bonds is encour- 
aging. With further development, this computationally 
tractable method holds great promise for complex or- 
ganic and biological systems where this manner of weak 
bonding plays a key role. 
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